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High  Order  Integration  of  Smooth  Dynamical  Systems  : 
Theory  and  Numerical  Experiments1 


By:  Mark  Austin 


Department  of  Civil  Engineering  and  Systems  Research  Center, 
University  of  Maryland, 

College  Park,  MD  20742. 


SUMMARY 

This  paper  describes  a  new  class  of  algorithms  for  integrating  linear  second  order  equations,  and  those 
containing  smooth  nonlinearities.  The  algorithms  are  based  on  a  combination  of  ideas  from  standard  New- 
mark  integration  methods,  and  extrapolation  techniques.  For  the  algorithm  to  work,  the  underlying  Newmark 
method  must  be  stable,  second  order  accurate,  and  produce  asymptotic  error  expansions  for  response  quan¬ 
tities  containing  only  even  ordered  terms.  It  is  proved  that  setting  the  Newmark  parameter  7  to  1/2  gives  a 
desirable  asymptotic  expansion,  irrespective  of  the  setting  for  (3.  Numerical  experiments  are  conducted  for 
two  linear  and  two  nonlinear  applications. 
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1  INTRODUCTION 


Engineers  need  fast  and  accurate  computational  methods  for  the  dynamical  design, 
analysis,  and  control  of  mechanical  and  structural  systems.  In  the  design  and  performance 
evaluation  of  dualspin  and  flexible  satellites,  for  example,  very  accurate  computational 
methods  are  needed  for  the  long  term  prediction  of  position  and  attitude.  Similarly,  fast 
computational  methods  are  needed  to  simulate  and  control  flexible  robotic  manipulators 
and  free  flying  robots.  Since  closed  form  analytic  solutions  are  unavailable  except  for 
a  few  specific  cases,  discrete  numerical  integrators  must  be  relied  on  to  predict  system 
behavior.  Early  numerical  experiments,  including  ours  [4,6],  exposed  weaknesses  in  the 
use  of  traditional  numerical  integrators.  The  use  of  standard  off-the-shelf  Runga-Kutta 
algorithms  to  integrate  the  dynamics  of  rigid  body  satellites  spinning  freely  in  space  nearly 
always  results  in  a  steady  accumulation  of  energy.  Thus,  they  are  a  poor  guide  to  the 
reliable  prediction  of  long  term  dynamical  phenomena,  and  highlight  the  need  to  design 
new  numerical  integrators  to  calculate  time-dependent  system  responses. 

To  provide  guidance  in  the  formulation  of  these  computational  models,  researchers 
are  attempting  to  identify  and  classify  the  natural  algebraic  and  geometric  structures  that 
these  dynamical  systems  possess  [15,17].  If  energy  and/or  norms  of  angular  momentum  are 
are  conserved  in  the  continuous  system  response,  then  it  is  very  desirable  the  same  enti¬ 
ties  remain  invariant  in  the  discrete  approximation.  In  the  case  of  canonical  Hamiltonian 
systems,  Feng  Kang  [14]  has  shown  that  it  is  possible  to  construct  symplectic  schemes  of 
arbitrary  order  accuracy,  which  are  A-stable  and  conserve  all  invariants  up  to  second  order. 
Some  recent  algorithms  [1,22]  for  noncanonical  Hamiltonian  systems  have  been  shown  to 
conserve  energy  exactly,  but  are  only  second  order  exact  in  their  response  estimates.  As  a 
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result,  when  large  time-steps  are  used  in  simulation  studies  of  long-term  dynamical  behav¬ 
ior,  displacements  can  be  completely  out  of  phase  after  only  several  hundreds  timesteps; 
examples  may  be  found  in  Austin,  Krishnaprasad,  and  Wang  [1],  and  Hoff  and  Taylor  [11]. 

Experience  indicates  that  while  numerical  algorithms  may  be  designed  to  conserve 
some  integrals  of  motion,  they  cannot  be  expected  (in  general)  to  conserve  all  integrals 
of  motion.  Marsden  et  al.  [17]  suggest  that  numerical  algorithms  be  crafted  to  conserve 
exactly  some  integrals  of  motion  important  for  design,  with  other  integrals  possibly  being 
conserved.  This  approach  has  several  problems.  First,  it  is  not  hard  to  find  engineering 
applications  where  the  equations  of  motion  are  so  nonlinear,  it  is  unlikely  that  a  discrete 
approximation  will  be  found  to  satisfy  even  one  invariant  of  motion.  And  when  the  system 
is  relatively  uncomplicated,  and  such  a  discrete  approximation  can  be  found,  performance 
may  still  be  less  than  satisfactory  (as  already  noted  above).  Moreover,  it  may  be  argued 
that  even  if  a  discrete  numerical  approximation  theoretically  conserves  selected  invariants 
exactly,  its  implementation  will  at  most  conserve  the  same  quantities  to  machine  precision; 
this  will  be  approximately  16  decimal  places  accuracy  for  double  precision  calculations  on 
standard  engineering  workstations. 

A  second  strategy  for  developing  numerical  algorithms,  and  the  approach  that  is  fol¬ 
lowed  in  this  paper,  is  to  find  efficient  ways  of  systematically  computing  displacements  and 
velocities  to  an  arbitrarily  high  order  of  accuracy.  This  approach  offers  several  computa¬ 
tional  benefits.  First,  it  is  easy  to  show  that  all  system  invariants,  which  Eire  simple 
polynomial  functions  of  displacement  and  velocity,  will  be  preserved  to  the  same  order  of 
accuracy  as  the  displacements  and  velocities.  Second,  the  equations  of  motion  may  be  nat¬ 
urally  partitioned  for  concurrent  computations  in  a  network  of  engineering  workstations; 
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for  details,  the  interested  reader  is  referred  to  Austin  and  Voon  [2]. 

The  purpose  of  this  paper  is  to  present  the  formulation  of  a  new  algorithm  that  uses 
ideas  from  Newmark  integration  and  extrapolation  techniques  to  accurately  compute  long 
term  dynamics.  The  scope  of  applications  is  limited  to  two  classes  of  finite  dimensional 
systems,  namely: 

[a]  Mass-Spring  Systems  :  Linear  undamped  (and  damped)  single  (and  multiple)  de¬ 
gree  of  freedom  systems  moving  under  free  and  forced  vibrations. 

[b]  Connected  rigid  body  assemblies  :  Applications  include  the  attitude  dynamics 
of  rigid  body  and  dualspin  spacecraft  subject  to  zero  external  forces,  dynamics  of 
floating  rigid  body  components  connected  by  hinges,  and  the  attitude  control  of  free- 
flying  robots.  For  discussions,  see  the  work  of  Sreeneth  [23,24],  Yang  [26],  and  Byrne 
[4],  Also,  see  Posbergh  [19]  for  the  formulation  of  a  geometrically  exact  rod  model  as 
a  series  of  springs  and  masses. 

Numerical  experiments  are  conducted  to  assess  the  effectiveness,  and  demonstrate  the  lim¬ 
itations  of  proposed  algorithm. 

2  EQUATIONS  OF  MOTION 

We  assume  that  the  dynamical  behavior  of  each  application  may  be  written  as  a  family 
of  n  2nd  order  differential  equations: 

M  [x(f),  x(f)]  x(t)  +  F  [x(t),  x(t)}  —  P  ( t )  (1) 

p 

with  initial  conditions  x(0)  and  x(0).  Here  x(t)  —  [xj(t),  X2(t)  ...  crn(t)]  is  a  (n  x  1)  vector 
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of  system  displacements,  M  [x(i),  x(f)]  a  (n  x  n)  mass-type  matrix,  F  [i(t),x(f)]  a  (n  x  1) 
general  force  vector,  and  P(t)  a  (n  x  1)  vector  of  external  loads  applied  at  the  nodal  degrees 
of  freedom.  An  alternative,  but  equivalent  form  of  (1)  is  obtained  by  setting  w(t)  =  i(i), 
and  rewriting  the  equations  of  equilibrium  as  2 n  first  order  equations,  i.e.: 

y(t)  = 

In  the  case  of  linear  second  order  dynamical  systems,  for  example,  often  F(w(t),x(t ))  = 
[c]rw(t)  +  [ k]x(t )  -  here  [c]  and  [k]  are  (n  x  n)  damping  and  stiffness  matrices  -  and 
M  [x(t),  x(t)]  is  a  constant  mass  matrix  [m].  Equation  (2)  may  be  written: 


x(t) 

w(t) 


=  )] 


w 


(*) 


M  1[w(t),x(t)}  ■  [ Pit )  -  F(w(t),x(t ))] 


(2) 


[v(01  = 


i(t)' 

>(<). 

o 

-i—i  i 


i  1 1 


x(t) 

w(t) 


+ 


0 

[m]-1[P(f)] 


(3) 


Equation  (3)  underscores  the  fact  that  the  complete  dynamics  are  characterized  by  time 
variations  in  displacements  and  velocities  alone.  Accelerations  are  recovered  via  back  sub¬ 
stitution  into  (1).  From  a  computational  standpoint,  it  is  desirable  to  work  with  (1) 
because  it  has  only  n  equations,  versus  2n  for  equations  (2)  and  (3). 


3  ONESTEP  INTEGRATION  SCHEMES 

Suppose  that  solutions  to  (l)-(3)  are  required  for  time  te  [tg,  f9+At]-  Irrespective  of 
whether  the  underlying  dynamics  are  described  as  families  of  first  order  or  second  order 
equations,  numerical  solutions  require  accurate  approximations  to: 


J/fc(fg+At)  — 


/fc(r>  y(r ))dr 


(4) 
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where  /fc(t,y(t))  is  the  kth  component  of  the  equation  of  motion.  In  problem  formulations 
based  directly  on  d’Alembert’s  principle,  [j/fc(t),  fk{t,  y(t))\  pairs  exist  for  displacement  and 
velocity,  and  velocity  and  acceleration,  respectively.  Similarly,  in  problem  formulations 
based  on  Lagrange’s  equations  and/or  Hamilton’s  equations,  [y*(i),  fk(t,  y(f))]  pairs  exist 
for  the  generalized  co-ordinates  and  momenta. 

A  wide  variety  of  onestep  numerical  integration  schemes  approximate  (4)  by  dividing 
the  time  interval  [tq,  t9+A<]  into  N  equal  intervals  of  distance  h  =  tq+k  —  tq,  and  sequentially 
applying  the  explicit  onestep  formula 

v(t,+h)  =  y(t,)  +  h  ■  *(*„»(<„),  h)  +  0(V+1).  (5) 

Here  $  :  Rx  H2n  x  1R  — >  1R2”  is  called  the  increment  function  of  the  method,  and  y(tq ) 
is  the  exact  solution  at  time  tq.  For  numerical  methods  that  have  implicit  increment 
functions,  we  write 

V(tq+h)  =  y(tq)  +  h-V  (tq,tq+h,y(tq),y(tq+h),h)  +  0(hp+1)  (6) 

where  ^  :  ]R  x  1R  x  IR2"  x  ]R2n  x  1R  — >  ]R2n.  It  is  important  to  note  that  when  the  onestep 
truncation  error  is  0(hp+1 ),  the  global  error  is  a  combination  of  local  plus  transported 
errors,  and  is  0(hp).  A  detailed  proof  is  given  in  Chapter  3  of  Hairer  et  al.  [9].  Moreover, 
multistep  methods  may  be  interpreted  in  the  framework  of  (4)  by  simply  rewriting  the 
numerical  schemes  as  onestep  methods  in  a  higher  dimensional  space  [9]. 
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4  BACKGROUND  TO  EXTRAPOLATION  METHODS 


Perhaps  the  most  straight  forward  way  of  increasing  numerical  accuracy  is  to  divide 
the  interval  [tg,ig+At]  into  smaller  increments  of  h.  While  this  strategy  may  be  acceptable 
for  short  periods  of  simulation,  it  places  the  designer  at  odds  with  the  need  to  compute  the 
long-term  time  evolution  of  dynamical  systems  without  resorting  to  excessive  computational 
effort. 

The  algorithm  developed  in  this  paper  hinges  on  the  existence  of  an  asymptotic  expan¬ 
sion  for  the  global  error  in  the  numerical  method.  Although  the  benefits  of  these  expansions 
were  known  to  Richardson  [20,21],  the  formal  theory  for  their  existence  is  due  to  Gragg  [8] 
and  Stetter  [25].  Gragg’s  ideas  are  most  succinctly  stated: 

(a)  Theorem  :  (due  to  Gragg)  Suppose  that  a  given  method  with  sufficiently  dif¬ 
ferentiable  (smooth)  increment  function  $(•)  satisfies  the  consistency  condition 
$(tg,  y(tq),  0)  =  /(fg,y(fg))  and  possesses  an  asymptotic  expansion: 

y(tg+fr)  =y(tq)  T  h  •  J/(tg), 

dp+1(t)hp+1  4-  •••  +  dp+k(t)hp+k  +  0(hp+k+1 )  (7) 

for  the  local  error.  If  N  steps  of  integration  are  computed  to  cover  the  time  interval 
At,  (i.e.  At  =  tq+At  ~  tq  =  N  ■  h),  then  the  global  error  has  an  asymptotic 
expansion  of  the  form: 

u(tg+At)  -  y(tg+At)  =  ep(t)  ■  A tp  4 - 1-  ep+k(t)  ■  A tp+k  4-  Eh{t)Atp  f  k+1  (8) 

where  ti(fg+At)  is  numerical  solution  at  time  At,  and  ep(t)  are  solutions  to  the 
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inhomogeneous  equations: 


e p (<)  =  -^(t,y(t))ep(t)  —  dp+1(t),  such  that  ep(tg)  =  0  (9) 

and  Eh(t)  is  bounded  for  tq  <t  <  tq+&t  and  0  <  h  <  h0.  /(•)  and  $(•)  are  as  defined 
in  equations  (2)  and  (5),  respectively. 

An  abbreviated  version  of  the  proof  to  Gragg’s  theorem  may  be  found  in  Hairier  amd  Lubich 
[11],  It  is  important  to  note  that: 

[a]  The  existence  of  equations  (7)-(9)  requires  /(•)  be  sufficiently  differentiable  over  the 
complete  time  interval  of  interest;  hence  use  of  the  word  smooth  in  the  title  of  the 
paper.  This  condition  automatically  precludes  the  use  of  extrapolation  methods  for 
applications  containing  sharp  material  (or  geometric)  discontinuities. 

[b]  As  stated,  equations  (7)-(9)  are  for  an  explicit  increment  function.  Subsequent  work 
by  Setter  [25]  has  established  that  asymptotic  expansions  of  the  form  (8)  also  exist  for 
implicit  onestep  methods. 

[c]  Extrapolation  techniques  require  the  underlying  numerical  method  to  be  stable.  Al¬ 
though  explicit  numerical  integrators,  such  as  Euler,  are  much  less  computational  work 
that  implicit  schemes,  they  are  notorious  for  being  only  conditionally  stable  [7,9].  In 
this  study,  base  steps  are  computed  with  implicit  increment  functions  that  are  A-stable. 

[d]  Preservation  of  the  asymptotic  expansion  requires  that  the  equations  (7)  be  solved 
exactly  at  each  timestep.  However,  it  is  demonstrated  in  the  following  sections  that 
even  when  nonlinearities  force  the  equations  to  be  solved  iteratively,  very  favorable 
increases  in  numerical  accuracy  are  still  possible. 
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The  most  important  extrapolation  methods  are  those  that  remain  invariant  under  the 


switching  of  parameters  tq  y(tq)  y(tq+h),  and  integration  order  h  —h.  For 

example,  if: 

y{iq+h)  =  y(iq )  +  h  ■  '&(tg,tq+h,y(tq),y(tq+h),h)  +  0(h3)  (10) 

is  equivalent  to 

V(iq)  =  V^g+h)  -  h  ■  ty(tq+h,tq,y(tq+h),y(tq),~h)  +  0(h 3)  (11) 

then  the  method  is  said  to  be  symmetric.  Instead  of  having  an  asymptotic  expansion  of 
the  form: 

U(tq+At)  =y(tq+At)  +  A2  At2  +  -43  At3  +  -4  4  At4  + 

4.5  At5  +  Aa  At6  +  -47  At7  H -  (12) 

where  the  coefficients  4.2,  4.3  •••  Aq  •••  are  solutions  to  (9),  numerical  integrators  with 
symmetric  increment  functions  have  asymptotic  expansions  of  the  form: 

U{tq+At )  =y(tq+At)  +  42  At2  +  -44  At4  +  46  At6  +  •  •  •  (13) 

Examples  of  implicit  numerical  methods  having  symmetric  increment  functions  are  the 
midpoint  rule 

y(t,+k)  =  y«,)  +  h  •  /(*->  *--«+»-,  (14) 

and  the  trapezoidal  method 

y(tq+h)  =  y(tg)  +  |  •  [/(t?,y(tg))  +  /(*9+A,y(*g+&))]-  (15) 
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Now  let  u™  (y(t9),  A t)  be  the  computed  numerical  solution  from  rik  steps  of  length  A f/n* 
using  a  numerical  method  of  mth  order  accuracy.  Furthermore,  assume  that  the  computa¬ 
tion  starts  from  exact  solution  y(tq).  Richardson  observed  [20,21]  that  improved  numerical 
approximations  could  be  obtained  by  solving  the  same  problem  over  a  prescribed  base 
step  A t  several  times  with  successively  smaller  internal  timestep  lengths.  If  the  under¬ 
lying  numerical  methods  is  2nd  order  accurate,  then  this  gives  a  sequence  of  numerical 
approximations: 

u2m  (y(*»),  A*),u£2  (y(tq),At),---,u2nk  (y(t,),At) 

for  the  set  of  increasing  integers  n\ ,  ri2  •  •  •  ru-.  Several  numerical  sequences  have  been  pro¬ 
posed,  including  that  of  Romberg  {1,  2, 4, 8, 16, 32,  •  •  •},  Bulirsch  [3]  {1,  2, 4, 6,  8, 12,  •  •  •}, 
and  the  double  harmonic  sequence  {2, 4, 6,  8, 10, 12, 14,  •  •  •}.  In  each  case,  linear  combina¬ 
tions  of  the  numerical  approximations  are  taken  to  eliminate  the  coefficients  A,-.  Succes¬ 
sively  refining  At  according  to  the  Romberg  sequence,  for  example,  gives: 

i=oo 

ui  (y(*j)>  At)  =  y(tq+&t)  +  Y  ^[At]2'  (16) 

i—  1 
1  =  00 

u2  (y(*g)»  Ai)  =  y(tq+At)  +  Y  A2t 

1=1 

Now  the  benefits  of  (13)  compared  to  (12)  are  evident  -  subtracting  (16)  from  4  times 
equation  (17)  gives: 

[4  •  y(tq+At)  -  y(tq+At)}  =  [4  ■  u|(-)  -  u .?(•)]  +  0  [At4] 

=  [4  —  1]  •  u4(-)  +  0  [At4]  (18) 

thereby  estimating  the  integral  with  two  additional  order  of  accuracy.  For  the  Romberg 
sequence  of  internal  timestep  lengths  it  is  relatively  straight  forward  to  show  that  successive 
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applications  of  the  formula: 


u 


(s+2) 

(2r) 


(0  = 


u(h') 


4(«/2)  _  1 


(19) 


will  systematically  eliminate  higher  order  coefficients  in  the  error  polynomial.  Here  r  is  the 
number  of  intervals  taken  to  integrate  across  the  A t  time-step,  and  s  the  order  of  accuracy 
for  the  numerical  estimate. 


5  NEWMARK  FAMILY  OF  INTEGRATION  METHODS 

Newmark  integration  methods  [12,18]  approximate  the  time  dependent  response  of 
linear  and  nonlinear  2nd-order  equations  by  insisting  that  equilibrium  be  satisfied  only  at 
a  discrete  number  of  points  (or  timesteps).  If  tq  and  tq+&t  are  successive  timesteps  in  the 
integration  procedure,  then  the  two  equations  of  equilibrium  that  must  be  satisfied  are: 

M  [i(tg),  *(*,)]  x(tq)  +  F  x(tq)}  =  P  (tq)  (20) 

and 

M  [x(tg r+A«)>  x(tq+&t )]  x(fg-j-At)  +  F  [x(tg+At),  ®(^g+A<)]  =  P  (^g+At)  •  (21) 


Now  let’s  assume  that  solutions  to  (20)  are  known,  and  (21)  needs  to  be  solved.  At  each 
timestep  there  are  3n  unknowns  corresponding  to  the  displacement,  velocity,  and  acceler¬ 
ation  of  each  component  of  x.  Since  we  only  have  n  equations,  the  natural  relationship 
existing  between  the  acceleration  and  velocity,  and  velocity  and  displacement  must  be 
enforced  to  reduce  the  number  of  unknowns  to  n.  That  is: 


®(^?+At)  —  X 


where  x(tq)  and  x(tq)  are  the  velocity  and  displacement  at  timestep  tq,  and  x(r)  is  an 
unknown  function. 

The  Newmark  family  of  integration  methods  assume  that:  (a)  acceleration  within 
the  timestep  behaves  in  a  prescribed  manner,  and  (b)  the  integral  of  acceleration  across 
the  timestep  may  be  expressed  as  a  linear  combination  of  accelerations  at  the  endpoints. 
Discrete  counterparts  to  (22)-(23)  for  the  update  velocity  and  displacement  are: 

x(tq+At)  =  x(tq)  +  At[(l  -  7 )x(tq)  +  7^(t9+A£)]  (24) 

At2 

x(tq+At )  =  x(tg)  +  A tx(tg)  +  — 2 - [(1  -  2 0)x(tq)  +  20x(tq+At)\  (25) 

with  the  parameters  7  and  0  determining  the  accuracy  and  stability  of  the  method  under 
consideration.  Equations  (24)  and  (25)  are  substituted  into  (21),  written  in  the  form 

g(-)  =  M[x(tq+At)]- x(tq+At)  -  P(<9+At)  =  0 


and  solved  by  iteration.  Finally,  x(tq+At)  is  back  substituted  into  (24)  and  (25)  for  the 
update  in  velocity  and  displacement. 


Remark  1  :  It  is  well  known  that  when  7  =  1/2  and  0  =  1/4,  acceleration  is  constant 
within  the  timestep  t  6  (f?,  tg+At],  and  equal  to  the  average  of  the  endpoint  accelerations, 
i.e: 


=  'Z(tq)  +  x(tq+At) 


=  x(tq)  + 


x(tq+At)  ~  x{tq) 

2 


=  x(tq)  + 


Ax(tg) 

2 


(26) 


In  such  cases,  approximations  to  the  velocity  and  displacement  will  be  linear  and  parabolic, 
respectively.  Moreover,  when  (20)  takes  the  form: 


13 


m]  x(tq)  +  [c]  x(tq )  +  [ k ]  x(tq )  =  P(fg) 


(27) 


equations  (21)  may  be  written 


[M]  A x(tq)  =  [AP] 


where 

and 


m  =  [m]  +  P-  [cl  +  [i] 

[AP]  =  P(t?+A«)  -  P(tq )  -  At  [c]  x(tq)  -  At  [fc] 


x(tq)  + 


(28) 


This  discrete  approximation  is  second  order  accurate  and  unconditionally  stable.  It  con¬ 
serves  energy  exactly  for  the  free  response  vibration  of  linear  undamped  SDOF  systems; 
see  Chapter  9,  page  512  of  Hughes  [12]  for  a  proof.  This  property  will  be  important  for  the 
algorithm  presented  in  the  following  sections. 
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6  NEWM ARK-EXTRAPOLATION  ALGORITHM 


The  algorithm  proposed  in  this  paper  draws  on  ideas  from  Newmark  integration  and 
Richardson’s  Extrapolation.  The  approach  is  motivated  by  a  desire  to  work  directly  with 
2nd  order  equations,  if  possible,  and  by  the  observation  that  improved  numerical  estimates 
of  response  will  result  when  equations  (22)  and  (23)  are  evaluated  with  increased  precision. 

Instead  of  simply  refining  the  step  length  of  the  numerical  method,  the  idea  is  to 
integrate  across  the  time  interval  t  €  several  times  with  successively  refined  step 

lengths.  In  each  case,  the  Newmark  method  is  used  as  the  base  integrator.  It  is  proposed 
that  with  a  judicious  choice  of  Newmark  parameters  7  and  (3  the  underlying  method  will 
be  second  order  accurate,  stable,  and  have  asymptotic  error  expansions  for  the  velocity 
and  displacement  components  containing  only  even  ordered  terms.  Thus,  extrapolation 
may  be  used  to  obtain  improved  estimates  of  displacement  and  velocity,  with  updates  in 
acceleration  being  computed  via  back  substitution  of  the  extrapolated  displacements  and 
velocities  into  the  equations  of  equilibrium. 

The  step  by  step  algorithm  for  P  levels  of  extrapolation  using  a  Newmark  base  method, 
parameter  settings  (3  and  7,  and  a  Romberg  sequence  of  refinement  is: 

[la]  Initialization  :  For  each  component  of  velocity  and  displacement,  dynamically  allo¬ 
cate  memory  for  a  (P  x  P)  extrapolation  table.  Lambert  [16]  reports  that  in  practice 
P  typically  falls  into  the  range  4-7. 

[lb]  Select  a  timestep  At.  In  some  cases  this  may  be  the  maximum  timestep  length  for 
algorithm  stability,  while  in  other  instances  it  may  be  a  suitably  small  At  needed  to 
draw  a  smooth  graph  of  response. 
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[2a]  Outer  Loop  :  For  Newmark  Integration,  q  =  1  to  nsteps 

[2b]  Set  tq  =  [q  —  1]  •  A t 

[3a]  Inner  Loop  1  :  For  i  =  1  to  P 

[3b]  Incremental  timestep  A U  = 

[4a]  Inner  Loop  2  :  For  k  =  1  to 

[4b]  Calculate  P(tk)  at  time  tk  =  tq  +  k  •  A ti 

[4c]  Solve  M[x(tk)]  •  x(tk)  =  P(tk) 

[4d]  Update:  =  x(tk-i)  +  Aij  •  [(1  -  j)x(tk-i)  +  /yx(tk)] 

x(tk)  =  x(tk-i)  +  At{  •  x(tk-i)  +  ■  [(1  -  2/3)x(tk-i)  +  2/3x(tk)\ 

[4e]  End  Inner  Loop  2 


[3d]  Put  displacement  x(tq+&t)  and  velocity  x(tq+ At)  components  in  position  uj(-)  of  Ex¬ 


trapolation  Table  [1]. 


Table  [1]  :  Extrapolation  Tableau  for  Romberg  Sequence 


[2c]  End  Inner  Loop  1  :  Note  :  Each  component  of  this  loop  may  be  computed  in 
parallel. 

[lc]  Use  equation  (19)  to  calculate  columns  2-P  of  each  extrapolation  table.  The  lower 
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most  right  entries  of  Table  [1],  u'[p  I^1^-),  are  taken  as  the  starting  displacement  and 
velocity  components  for  the  next  time-step. 

[ld]  Back  substitute  the  displacement  and  velocity  vectors  from  Step  [lc]  into  equation 
(21),  and  solve  for  acceleration  vector  £(t9+A<)- 

[le]  End  Outer  Loop 

The  Newmark  parameters  7  and  /9  are  selected  according  to  a  dual  set  of  criteria;  in 
addition  to  requiring  the  base  method  be  second  order  accurate  and  stable,  asymptotic 
error  expansions  for  updates  in  velocity  and  displacement  must  contain  only  even  ordered 
terms. 

6.1  Accuracy  and  Stability 

For  linear  dynamical  systems  of  the  type  mentioned  above,  the  Newmark  method  is 
second  order  accurate  if  and  only  if  7  =  1/2.  Setting  7  =  1/2  and  /?  =  1/4  results  a 
method  that  is  unconditionally  stable.  For  nonlinear  systems  y  =  f(t,  y),  which  are  twice 
continuously  differentiable  in  t  and  y,  the  trapezoidal  rule  is  stable,  convergent,  find  second 
order  accurate  [13].  The  other  common  choice  of  parameters  is  7  =  1/2  and  ft  ~  1/6;  it 
corresponds  to  a  linear  time  variation  in  acceleration,  quadratic  variation  in  velocity,  and 
cubic  variation  in  displacement.  The  price  for  increased  accuracy  (i.e.  a  smaller  coefficient 
for  the  local  truncation  error)  is  conditional  stability.  This  feature  is  of  little  benefit 
because  the  purpose  of  the  extrapolation  is  to  eliminate  the  truncation  errors  completely, 
irrespective  of  their  size. 
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6.2  Asymptotic  Error  Expansion 

The  goal  of  this  section  is  to  determine  which  settings  of  7  and  ft  result  in  asymptotic 
expansions  for  the  local  error  containing  only  even  ordered  powers  of  At.  First,  assume 
that  M  [•]  is  invertible.  Equations  (20)-(21)  may  be  written: 

x(tq)  =  M_i  x(tq)]  •  [P  (tq)  —  F  [x(tq),  x(tq)}}  (29) 

and 

x{tq+At)  =  M-1  [i{tq+At),x{tq+At)]  ■  [P  {tq+At)  ~  F  [i(tg+  At  ),  x(tq+  At)]]  (30) 

Substituting  equations  (29)-(30)  into  (24)-(25)  and  rearranging  terms  gives  a  pair  of  implicit 
equations  for  the  update  in  velocity  and  displacement.  They  are: 

X  {tq+  At )  =  ^{tq)  T  At  •  \F  i  {tq,tq^.Ati  x{tq)-i  x{tq),  x{tc jr-f  Ai  )  i  ^{tq+At)i  At )  (31 ) 

where  'Fi  (•)  =  (1  -  7)  •  M~l  [i(t3),  x(tq))  •  [P  (tq)  -  F  [x{tg),  x(<?)]]  + 

7  •  M-1  [x{tq+At),x{tq+At)}  '  [F  {tq+At)  ~  F  [i(ig+At),  x{tq+At )]] 

and 

x{tq+At  )  — -  x{ tq )  d-  At  •  *F  2  {tqt  tq+At  i  x{t  q) ,  x{t  q) ,  x{t  q-\-  At  )i  %{tq+  At  )i  At)  (3— ) 

where  *2(.)  =  i(ig)  +  ^  •  (1  -  2ft)  ■  M~l  \x{tq),  *(*,)]  •  [P  (tq)  -  F  [i{tq),  *(*,)]]  + 

•  2ft  •  M  1  [x{tq-\-At)i  X{tq+At  )]  '  [P  {tq+At)  F  [x^tq-ij-  At)  i  X{tq+At)}}  ■ 

Assume  that  'Fi(-)  and  \F2(-)  are  continuously  differentiable  with  respect  to  the  displace¬ 
ments  and  velocities  over  the  complete  time  interval  [t?,tg+A<]-  Equations  (31)  and  (32) 
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will  have  even  ordered  asymptotic  expansions  if  the  interchange  of  parameters  tg  <-»  tq+At > 


x{tq)  <->  :r(tg+At),  i(tg)  i(tg+At)  and  At  —At  leaves  the  equations  invariant.  In 
other  words,  equation  (31)  must  be  identical  to 

x{tq  )  =  &  {tq+At )  At  •  \&i  {tq+At  itq,  x{tq+At  )j  {tq+At  ),3;(tg),i(tg),,  At)  (33) 
where  $1  (•)  =  (1  -  7)  •  M_1  [i(tg+ A<),  z(/9+At)]  •  [-P(ig+At)  -  -F  [z(tg+A1 ),  ®(<?+At )]]  + 

7  •  ^_1  [*(*?)>  *(*g)]  •  [-P(<9)  -  F[i(tg),x(tg)]] 

Equating  coefficients  gives  1  —  7  =  7,  he.  7  —  1/2-  Similarly,  the  displacements  will 
have  even  ordered  asymptotic  expansions  if  the  interchange  of  parameters  tg  <-4  tg+A*, 
x(tg)  <-4  x{tq+At),  i{tq)  <-»  i(tg+At)  and  At  «-»  -At  in 

At2 

x(tg+At)  =  x(tg)  +  A ti(tg)  +  -y-[(l  -  2 0)x(tq)  +  2/3i(tg+At)]  (34) 

gives  a  numerical  scheme  equivalent  to 

At2 

x(tq)  =  x(tq+At )  -  Ati(tg+At)  +  -^-[( 1  -  2/3)x(tg+A<)  4-  2/?x(tg)]  (35) 

Recall  that  7  =  1/2.  Substituting 

i(tg+At)  =  x(tq)  +  ~-[x{tq)  +  x{tq+At)] 

into  (35)  and  rearranging  terms  gives 

X{tq+At)  —X(tq)  +  A  tx(tq)  + 

At2  At2 

— [i((,)  +  i((,+A,)] - 2~[(1  -  2/J)i(i,+A.)  +  2/3i(f,)].  (36) 

It  is  evident  that  the  coefficients  in  (34)  and  (36)  will  be  identical  for  any  value  of  (3.  ■ 
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Remark  2:  A  natural  question  to  ask  is  whether  or  not  “other  combinations  of  7  and  0 
are  admissible  when  the  system  dynamics  are  linear  ?”  Consider,  for  example,  the  linear 
SDOF  model 

x(tq)  +  2  £U!x(tq)  +  UJ2x(tq)  =  P(tq). 

For  the  purposes  of  conducting  the  analysis,  it  is  convenient  to  define  the  column  vector 
z(tq)  —  [x(fg),  x(tq)]T  and  write  the  discrete  update  for  Newmark  in  first  order  form 

A\  [At]  •  2(tg4-A«)  =  A2  [At]  •  z{tq )  +  L[tq,tq+At,  At]  (37) 


where  A\  [At]  = 


1  +  A  t20uj2  2  At2  £7 


U) 


Atju2  1  +  2eAtjuj 


A2  [At]  = 


1  -  At2/2(1  -  20)u>2  At(l  -  At(l  -  2 0)eu) 
— At(l  —  j)u)2  1  —  2At(l  —  7  )eu 


and  L  [tg,  tg+A«,  At]  = 


A(2/2.[(l-2/3)P((5)  +  WP(t,+A,)\ 
At  •  [(1  -t)P(t,)  +  7.P(f,+A1)] 


Recasting  equations  (37)  in  the  form  of  (6)  gives 


~(^g-|-Ai)  —  z{tq)  4"  At  •  (tg,  tg-|-AA  ^(tg),  At) 


(38) 


where  ^  1  ( • ) 


'VM 

At 


•[(A2[At] 


^4-1  [At])  •  2r(tg)  +  L[tq,  tg+A«j  At]] .  (39) 


In  is  important  to  note  that  in  (38)  the  notation  'I'i(')  is  used  for  a  (2n  x  1)  column  vector 
composed  of  equations  'J'i(-)  and  vE^-)  in  (31)  and  (32),  respectively.  The  interchange  of 
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parameters  tg  <-*•  tg-f  At  <->■  —  At,  and  ;z(tg)  <->  z(tg-j-At)  will  leave  the  numerical  method 
unchanged  if  and  only  if 

A"1  [At]  •  A2[At]  =  A-'i-At]  •  [-At]  (40) 

and  L[tg,tg+A«,  At]  +  A\  [At]  •  A?1  [-At]  •  X[t?+At,  tg,  -  At]  =  0.  (41) 

Plugging  equations  (40)-(41)  into  Mathematica  [27]  and  solving  for  solutions  gives  7  =  1/2 
and  any  value  (3. 
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7  NUMERICAL  EXPERIMENTS 


This  section  presents  the  results  of  four  numerical  experiments;  three  examples  and 
one  counter  example.  Unless  otherwise  stated,  simulations  are  conducted  using  7  =  1/2 
and  (3  —  1/4.  The  extrapolation  computations  consists  of  four  levels  of  stepwise  refinement 
matching  the  Romberg  sequence,  i.e.  At,  At/2,  At/4,  and  At/8,  i.e.: 


I 


zs.t 


^t/2 


X 


zp.  t/4 


\ 


^t/a 

_ i\ 


i 


1  step 


2  sraps 

4  steps 


a  steps 


They  are  conducted  without  stepsize  and  error  control  in  the  sense  of  Deuflhard  [7].  In  order 
to  benchmark  the  accuracy  of  the  proposed  Newmark-Extrapolation  algorithm  against 
traditional  approaches  (for  equivalent  computational  work  per  unit  increment  of  time),  the 

timestep  for  the  standard  Newmark  integration  procedure  is  divided  by  15  -  i.e.  step  length 

At 
15  ’ 

The  accumulation  of  differences  between  analytic  and  numerical  response  qiiantities  is 
captured  by  defining  the  area  of  component  error  as: 

i=N 

Error  Area  =  stepsize  •  jAnalytic(f)  —  Numerical^)!;  (42) 

i= l 

where  N  equals  the  total  number  of  integration  timesteps  over  the  simulation  time  inter¬ 
val.  The  parameters  Analytical{t )  and  N umerical{t)  are  analytic  and  numerical  solutions 
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Component  x  :  Analytic  x(0.03)  =  -0.47884882915567743984 

Steps 

Stepsize 

0  [At2] 

0  [At4] 

0  [At6] 

00  ' 
-to 

□ 

1 

At 

-0.47827819849 

2 

At 

2 

-0.47870594155 

-0.47884852258 

4 

At 

4 

-0.47881309285 

-0.47884880995 

-0.47884882911 

8 

At 

8 

-0.47883989418 

-0.47884882795 

-0.47884882915 

-0.478848829155675 

Component  x  :  Analytic  i(0.03)  =  0.99280863585386625224 

Steps 

Stepsize 

o  [At2] 

0  [At4] 

0  [Ai«] 

□  [At8] 

1 

At 

0.99282582702 

2 

At 

2 

0.99281294252 

0.9928086477 

4 

At 

4 

0.99280971308 

0.9928086366 

0.99280863586 

8 

At 

8 

0.99280890519 

0.992808635853866 

Table  [2]  :  Simulation  components  after  1  Time  Step 


to  components  of  response  (acceleration,  velocity,  and  displacement),  and  system  invari¬ 
ants  (total  angular  momentum,  energy  and  so  on).  All  calculations  were  done  in  double 
precision  arithmetic  on  a  SUN  workstation. 

7.1  Example  1  :  Undamped  SDOF  Oscillator 

The  free  vibration  response  of  an  undamped  SDOF  oscillator  is  studied  as  a  means  of 
gaining  insight  into  the  behavior  and  properties  of  the  Newmark-extrapolation  algorithm. 
If  the  mass  m  =  1  and  stiffness  k  =  16,  with  initial  conditions  x(0)  =  1  and  x(0)  =  0,  then 
the  analytic  solution  is: 


x(t)  —  cos(4t) 


(43) 
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Propagation  of  Energy  Errors  | 

Steps 

Stepsize 

0  [At2] 

0  [At4] 

0  [At6] 

CO 

-V-3 

<L 

□ 

1 

At 

0.000000 

2 

A  t 

2 

1.77e-15 

4.12e-08 

4 

At 

4 

1.77e-15 

2.58e-09 

l.lle-11 

8 

At 

8 

-3.55e-15 

1.61e-10 

1.66e-13 

-6.21e-15 

Table  [3]  :  Energy  Errors  for  Newmark-Extrapolation  at  Time  Step  1 


The  simulation  timestep  in  the  Newmark-Extrapolation  algorithm  was  set  at  At  =  0.03 
seconds.  For  equivalent  computational  work  per  unit  time,  the  simulations  were  repeated 
for  the  standard  Newmark  algorithm  using  At  =  0.002  seconds  and  15  times  the  number 
of  time  steps. 

Performance  of  Numerical  Algorithm  at  Step  1 

Tables  [2]  summarizes  the  component  response  of  the  Newmark-Extrapolation  algo¬ 
rithm  after  only  1  timestep;  i.e.  t  =  0.03  seconds.  The  entries  of  column  3  correspond  to 
velocity  and  displacement  responses  computed  with  the  standard  Newmark  Method.  Once 
column  3  is  filled  in,  equation  (19)  is  used  for  the  calculation  of  columns  4-6,  with  elements 
Ug(-)  being  taken  as  the  starting  displacement  and  velocity  at  timestep  2.  A  key  point  to 
note  from  Table  [2]  is  that  convergence  of  the  numerical  solution  to  the  analytic  solution 
is  much  faster  due  to  extrapolation  (across  the  rows  of  the  table)  than  by  reduction  of  the 
step  size  (moving  down  the  columns). 

Physical  considerations  dictate  that  energy  will  be  conserved  in  (43).  Table  [3]  tracks 
the  errors  in  the  system  energy  at  each  stage  of  the  extrapolation  process  for  time  step  1. 
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As  indicated  in  the  opening  sections,  the  standard  Newmark  method  with  7  =  1/2  and  j3  = 
1/4  theoretically  conserves  energy  exactly  for  this  application,  irrespective  of  the  timestep 
length.  In  practice,  however,  minor  errors  are  introduced  due  to  numerical  roundoff  (see 
column  3  of  Table  [3]).  Column  4  of  Table  [2]  contains  displacement  and  velocity  estimates 
that  are  4th  order  accurate.  The  corresponding  column  in  Table  [3]  contains  large  errors 
in  the  energy  of  the  numerical  response;  this  observation  is  consistent  with  Dahlquist’s 
stability  criterion  for  this  integration  linear  systems  using  a  multi-step  method  [5].  As  the 
extrapolation  process  continues,  progressively  higher  order  estimates  of  response  quantities 
are  obtained,  and  the  energy  error  systematically  approaches  zero.  Indeed,  the  final  error 
in  energy  is  of  the  same  order  as  would  occur  due  to  numerical  roundoff. 

Long  Term  performance  of  Numerical  Algorithm 

The  numerical  response  of  the  Newmark-Extrapolation  algorithm  was  calculated  for 
3000  seconds,  using  to  100,000  time  steps  of  At  =  0.03  seconds.  This  corresponds  to 
approximately  1910  full  cycles  of  (43). 

The  first  block  of  simulations  corresponds  to  parameter  settings  7  =  1/2  and  j3  =  1/4. 
The  Newmark-Extrapolation  algorithm  (see  column  3  of  Table  [4])  tracks  the  time  varia¬ 
tion  in  acceleration,  velocity  and  displacement  quantities  much  closer  than  the  standard 
Newmark  method.  The  error  area  in  energy  for  the  standard  Newmark  method  should  the¬ 
oretically  be  zero.  However,  an  accumulation  in  roundoff  errors  gives  an  absolute  integral 
of  energy  errors  4.712  x  10-11.  The  corresponding  error  area  in  energy  given  by  Newmark- 
Extrapolation  is  of  the  same  order  as  the  displacements  and  velocity,  as  expected. 

A  second  set  of  simulations  were  conducted  using  7  =  1/2  and  (5  =  1/6.  The  objective 
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Error  Component 

Standard  Newmark 

Newmark  Extrapolation 

Displacement 

4.074  x  10° 

5.304  x  10-8 

Velocity 

1.629  x  101 

2.121  x  10-7 

Acceleration 

6.518  x  101 

8.487  x  10~7 

Energy 

4.712  x  10-11 

9.514  x  10~8 

Table  [4]  :  Error  Area  after  3000  Seconds  :  gamma  =  1/2  :  beta  ==  1/4. 


Error  Component 

Standard  Newmark 

Newmark  Extrapolation 

Displacement 

2.037  x  10° 

8.474  x  10~9 

Velocity 

8.148  x  10° 

3.389  x  10~8 

Acceleration 

3.259  x  101 

1.356  x  10"7 

Energy 

4.267  x  10“3 

1.186  x  10~8 

Table  [5]  :  Error  Area  after  3000  Seconds  :  gamma  =  1/2  :  beta  ==  1/6. 


was  to  verify  that  extrapolation  would  work  for  values  of  (3  other  than  1/4,  and  to  see  how 
much  a  linear  variation  in  acceleration  within  the  time  step  -  versus  constant  acceleration 
-  affects  numerical  accuracy.  In  the  analysis  of  a  linear  undamped  free  vibrating  SDOF 
systems,  it  can  be  shown  -  see  Chapter  9  of  Hughes  [12]  -  that  this  Newmark  method  does 
not  conserve  energy  exactly,  and  is  stable  only  when  u>  ■  A t  <  3.464.  For  this  application, 
the  prodiict  io  •  At  =  0.12,  so  stability  is  not  a  problem. 

The  results  of  Tables  [4]  and  [5]  indicate  that  a  linear  approximation  in  accelera¬ 
tion  improves  the  numerical  accuracy  in  displacements,  velocities,  and  accelerations  in 
both  the  standard  Newmark  method  and  the  proposed  Newmark-Extrapolation  algorithm. 
For  example,  estimates  of  displacement  from  Newmark-Extrapolation,  have  error  areas 
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53.04  x  10-9  and  8.474  X  10-9  for  the  constant  and  linear  acceleration  approximations, 
respectively.  This  improvement  in  numerical  accuracy  is  minor,  compared  to  that  obtained 
via  extrapolation. 

Comparison  to  Mid-point  Rule 

In  the  numerical  solution  of  complicated  dynamical  systems,  often  the  equations  of 
motion  are  written 


g(-)  =  M  [x(tfl+At)]  •  x(tq+&t)  ~  P(tg+At)  =  0 


and  solved  by  iteration  at  each  timestep.  One  problem  with  this  approach  is  the  potential 
destruction  of  the  asymptotic  error  expansion  if  the  iterative  solution  to  the  equations 
is  not  sufficiently  accurate;  that  is,  compared  to  the  so-called  exact  method  using  LU 
decomposition.  To  see  if  this  was  likely  to  happen,  the  SDOF  dynamics  were  written 

(44) 

Here  x(t)  and  p(t)  are  the  displacement  and  momentum  -  p(t)  =  m  •  x(t )  -  of  the  mass- 
spring  system.  The  discrete  approximation  to  (44)  using  the  midpoint  rule  is 


p(t) 

‘0  -1‘ 

r i  oi 

m 

.*(*). 

1  0  _ 

0  k. 

x(t)_ 

p(tq+6.t)-p(t7) 

At 

Z(t<,4-At  )~x(tg) 
At 


'0  -1' 

r i  oi 

771 

1  0 

0  k 

P(<<j  +  At  l  +  Pij^) 
2 

x(tq+At)  +  x(tq) 
2 


(45) 


The  alert  reader  will  notice  that  energy  is  conserved  when  equations  (45)  are  solved  exactly. 
Furthermore,  in  the  solution  of  linear  dynamical  systems,  the  midpoint  and  trapezoidal 
rules  are  equivalent  simply  because: 

'*(*,)  +  Z(i,+A,p  _  [/(*(,,))  +  /(l((,+  A,))] 


At-/ 


V 
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Error  Component 

Newmark- Extrapolation 

Midpoint-Extrapolation 

Displacement 

5.304  x  10“8 

5.222  x  10~8 

Velocity 

2.121  x  10“7 

2.089  x  10~7 

Energy 

9.514  x  10~8 

9.498  x  10-8 

Table  [6]  :  Error  Area  after  3000  Seconds 


Numerical  solutions  to  (45)  were  calculated  for  100,000  timesteps  of  At  —  0.03  sec¬ 
onds,  the  equations  being  solved  at  each  timestep  using  the  iterative  procedure  described 
in  Section  7.3.  Table  [6]  shows  the  error  areas  for  the  Newmark- Extrapolation  and  the 
midpoint  formulation  are  similar,  a  good  indication  that  the  iterative  equation  solver  is 
working  well. 

7.2  Example  2  :  Damped  SDOF  subject  to  External  Loading 

Consider  the  damped  SDOF  subject  to  external  loads: 


x(t)  +  4  x(t)  +  13x(t)  = 


-2 t 


sin[3t] 


(46) 


If  x(0)  =  1  and  x(0)  =  —2,  then  Laplace  transforms  gives  an  analytic  solution: 

x(t)  =  e~2t  •  cos[3t]  +  [— — —  I  [sin(3t)  —  3fcos(3i)]. 
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(47) 


The  damped  natural  period  of  this  system  is  approximately  2.3  seconds. 

Simulation  Results 

200  timesteps  of  simulation  were  computed  using  Newmark  parameters  7  =  1/2,  (3 
=  1/4,  and  timestep  length  At  —  0.03  seconds.  This  corresponds  to  approximately  3  full 
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cycles  of  motion.  Table  [6]  shows  that  the  standard  Newmark  algorithm  produces  error 
areas  of  the  order  10~6.  When  the  Newmark-Extrapolation  algorithm  is  used  to  compute 
the  response,  the  same  error  areas  are  of  the  order  10~15. 


Error  Component 

Standard  Newmark 

Newmark  Extrapolation 

Displacement 

2.32  x  10“6 

8.067  x  10~16 

Velocity 

9.30  x  10-6 

3.223  x  10-16 

Acceleration 

3.36  x  10"6 

1.235  x  10-15 

Table  [7]  :  Error  Area  after  200  Time  Steps 


7.3  Example  3  :  Lagrangian  Formulation  for  2-Body  Problem 

The  third  application  is  a  study  of  the  dynamical  behavior  of  two  planar  rigid  bodies 
connected  by  a  frictionless  revolute  joint  in  a  zero-gravity  environment.  Let  the  bodies 
have  masses  mi  and  m2,  and  moments  of  inertia  about  their  centers  of  mass  I\  and  I2. 
Figure  [1]  shows  that  the  centers  of  mass  for  each  body  are  located  at  distances  d\  and  d2 
from  the  revolute  joint.  Moreover,  the  orientation  of  each  body  is  described  by  the  angle  - 
8\{t)  and  #2 (t)  -  it  makes  to  the  x-axis,  measured  in  an  anti-clockwise  direction. 


In  the  absence  of  external  torque  the  system  is  conservative.  The  total  angular  mo¬ 
mentum  of  the  system  is 


M  =  fl  1 


Ji  e  •  d\d2  ■  cos(62(t)  —  #1  (t)) 

e  •  d\d2  ■  cos(d2(t)  —  I2 


ex{t) 

hit) 


.  (48) 


and  the  Lagrangian  of  the  system 
_  1 


IxO\(t)  -1-  I29l(t)  4-  ■  92(t)e  ■  dx d2  ■  cos(92(t)  —  #i(t)) 


(49) 


29 


Figure  [l]  Configuration  of  Planar  Two  Body  Dynamical  System 


In  equations  (48)  and  (49),  e  =  (mim2)/(mi  +  m 2)  is  a  reduced  mass,  and  I\  ==  I\  +  ed\ 
and  I2  =  I2  +  ec?2  are  augmented  inertias  of  the  bodies.  The  Euler-Lagrange  equation 
is  used  to  .derive  (see  Sreenath  et  al.  [23,24]  for  details)  the  equations  of  motion  for  this 
system.  They  are: 


edi d2  ■  cos(02 (t)  ~  (f ))  1 

"m 

_edid2  ■  cos (02 (0  —  #i(t)) 

- 1 

W*)\ 

edid2  ■  sin (92(t)  -  9i{t)) 


-922(t) 

m 


= 0. 


(50) 


Notice  that  (50)  is  linear  in  angular  accelerations,  and  falls  into  the  class  of  problems 
covered  in  (1).  If  angles  and  velocities  are  provided  as  initial  conditions,  then  the  starting 
angular  accelerations  may  be  recovered  directly  from  (50). 
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Computational  Implementation 

The  first  step  of  the  computational  procedure  is  to  write  (50)  in  component  form, 
namely: 

gi(di(t),92(tjS)  =  hOi (t)  +  e-  d\d2  ■  cos(92(t)  ~  Qi(t))92(t) 

—  B\{t)  ■  ed\d2  ■  sin(92(t)  —  =  0,  (51) 

92  (jii(t)i  92 (f)^  =  ^2^2  +  e  •  d\d2  •  cos( 92(t)  —  9i(t))d\(t) 

+  9\{t)  •  e  •  did2  ■  sm(82(t)  -  9-i(t))  =  0,  (52) 


where  the  subscripts  take  there  usual  meaning.  The  number  of  unknown  variables  at  each 
timestep  is  reduced  to  two  by  prescribing  constant  angular  acceleration  across  each  time 
interval  t  €  At]-  In  an  analogous  manner  to  equations  (22)-(23),  it  follows  that: 


.  .  rti+*‘  ..  .  At  r- 

0(t,+At)  =  6(t,)  +  J  0(T)dT  =  e(t,)  +  —  0(t,)  +  8(t,+ii) 


(53) 


and 


9(tq+At)  —  9q  + 


rb+At 


At 


0(r)dr  =  9(tq)  +  9(tq)At  q — —  0(tg)  +  6(tq+ At)  •  (54) 


Substituting  (53)-(54)  into  (51  )-(52)  gives  two  nonlinear  equations  where  the  only  unknowns 
are  #(i9+At)-  An  iterative  solution  is  computed  by  letting: 


«Sa,  =  +  Aflf. 


i>2+,U,  =  «2.,+A,  +  A9J. 


(55) 

(56) 


where  9^  q+At  is  the  kth  iterate  for  acceleration  component  i  at  time  here  we  have 

dropped  the  notation  0(fg+At)  only  to  make  the  equations  more  compact.  A  sequence  of 
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iterative  solutions  is  defined  by  taking  a  Taylor  series  expansions  about  g\  #2  ?+  At) 

and  <72(^1  q+Av  @2  q+At)  an-d  truncating  higher  order  terms.  This  gives: 


=^0tq+Ajiq+A  t) + 1  +  ^'2  + 

^1,9+At  C)^2,q+At 


O((A0*)2,(A02fc)2)  -  0, 


=S2(«f,!+A,.«2,,+A,)  +  T-A«f  +  dr^-A**  + 

+At 


72,g+A< 


O((A0*)2,(A0*)2)  =  0. 


Truncating  the  nonlinear  terms  of  A0  in  (57)-(58),  and  rewriting  in  matrix  form  gives 


#!  (01,9+ At’  @2, q+At) 
#2(01, 9+ At  5  0*2, 9+ At)  - 


(59) 


which  is  solved  for  the  incremental  update  in  components  (55)-(56).  A  concise  statement  of 
(59)  is  J.h  =  -g.  Iterations  continue  at  each  timestep  until:  (a)  a  preset  maximum  number 
of  iterations  is  reached,  or  (b)  all  changes  in  angular  acceleration  components  from  (55)- 
(56)  are  less  than  a  preset  error  value  times  the  magnitude  of  the  acceleration  imbalance  at 
the  beginning  of  the  iterations.  Moreover,  divergence  of  the  iterates  is  avoided  by  ensuring 

lb(01,t+At>02,t+At)ll2  is  less  than  II #(01,9+ At’  02,9+ At) II 2 '  When  this  test  fails’  h  is  divided 
by  powers  of  2  until  the  inequality  is  satisfied. 


Simulation  Results 

Consider  the  response  of  the  2-body  having  properties  m\  =  1,  m2  =  2,  d\  =  1.0,  c?2 
=  1.5,  Ii  =  1.0,  and  I2  =  3.0.  If  the  initial  displacement  and  velocity  vectors  are  [0, 1]  and 
[0.0,  5.0],  respectively,  then  back  substituting  into  equation  (50)  gives  (approximate)  start- 
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Error  Component 

Standard  Newmark 

Proposed  Algo 

Lagrangian 

4.421  X  1CT3 

9.309  x  10~10 

Momentum  Norm 

1.351  x  10“3 

2.293  x  10_1° 

Table  [8]  :  Error  Area  after  1000  Time  Steps 


ing  angular  acceleration  components  [13.133,  —1.5768].  The  norm  of  angular  momentum  is 
25.20151152934070  and  the  Lagrangian  56.25. 

The  response  of  the  2-body  was  calculated  for  1000  timesteps  at  At  =  0.03  seconds 
using  the  Newmark-Extrapolation  algorithm  (i.e.  t  €  [0,30]  seconds).  The  computation 
was  repeated  using  the  Standard  Newmark  algorithm  for  15000  timesteps  at  At  —  0.03/15 
=  0.002  seconds.  In  both  cases  this  interval  of  simulation  corresponds  to  approximately  10 
full  cycles  of  the  coupled  rigid  body  components. 

As  indicated  in  equations  (48)  and  (49),  theoretical  considerations  dictate  that  the 
body  momentum  and  energy  will  be  constant.  Consequently,  error  areas  for  momentum 
and  energy  were  calculated  for  both  simulation  cases,  and  are  shown  in  Table  [8].  They 
indicate  that  the  Newmark-Extrapolation  is  significantly  more  accurate  than  the  standard 
Newmark  algorithm. 

8.4  Counter  Example  4  :  Bilinear  Mass-Spring  System 

The  counter  example  demonstrates  how  the  extrapolation  process  will  fail  when  the 
underlying  equations  of  motion  are  not  sufficiently  differentiable,  as  required  by  Gragg’s 
theorem.  Consider,  for  example,  the  free  vibration  of  a  simple  mass-spring  SDOF  oscillator 
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having  mass  m,  and  a  bilinear  spring  restoring  force 


F(x) 


10m,  zone  1  if  |  x  |  <  1.0; 

10  *  (1  +  p(x  —  1)),  zone  2  for  |  x  |  >  1.0. 


where  x  is  the  displacement  of  the  mass  and  p  is  a  coefficient  of  strain  hardening.  If  the 
mass  has  displacement  a;(0)  and  velocity  i(0)  at  time  t  =  0,  then  the  analytic  solution  for 
free  vibration  is  a  piecewise  series  of  solutions 


x(t)  =  Acos(w(x)t)  +  Bsin(w(x)t) 

for  displacements  within  each  zone.  The  instantaneous  circular  frequency  of  the  system 
is  given  by  w(x)  =  y/{k(x)/m\,  where  k(x)  =  dF{x)/dx  is  the  tangent  stiffness.  The 
coefficients  A  and  B  depend  on  the  initial  conditions  as  the  system  passes  from  one  zone  to 
another.  More  important,  let  T(-)  be  the  increment  function  -  in  the  sense  of  equation  (6) 
-  for  the  numerical  method  when  the  trapezoidal  rule  is  used  to  approximate  the  update 
in  velocities  and  displacements.  There  does  not  exist  a  Lipschitz  constant  L  such  that 

\  gy.fefh  -  i  <  iin  -*2| 


for  all  x\  and  x 2  in  the  neighborhood  of  x  =  ±1.  From  a  physical  viewpoint,  this  translates 
to  a  discontinuity  in  the  third  derivative  in  displacement  at  the  boundary  between  zones 
1  and  2.  The  necessary  conditions  for  Gragg’s  theorem  are  not  satisfied,  and  as  a  result, 
the  extrapolation  should  fail  as  the  boundary  between  zones  1  and  2  is  crossed. 

Simulation  Results 

For  the  numerical  experiments,  mass  m  was  set  to  1,  strain  hardening  ratio  p  set  to  0.5, 
and  the  initial  displacement  and  velocity  at  t  =  0  set  to  2  and  0,  respectively.  The  system 
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Component 

Analytic 

Numerical 

Error 

Displacement 

1.079367284205931 

1.079367284205930 

0.000000000000001 

Velocity 

-4.835406755059004 

-4.835406755059001 

-0.000000000000003 

Acceleration 

-10.396836421029658 

-10.396836421029651 

-0.000000000000007 

Energy 

17.500000000000000 

17.499999999999972 

0.000000000000028 

Table  [9]  :  System  Components  at  end  of  Timestep  9 


Component 

Analytic 

Numerical 

Error 

Displacement 

0.877842947098626 

0.877954601994794 

-0.000111654896169 

Velocity 

-5.227371560478852 

-5.224170476570483 

-0.003201083908370 

Acceleration 

-8.778429470986259 

-8.779546019947945 

0.001116548961686 

Energy 

17.50000000000000 

17.499999999954522 

0.000000000045478 

Table  [10]  :  System  Components  at  end  of  Timestep  10 


has  constant  energy  =  17.5.  If  A t  =  0.04  seconds,  then  the  Newmark- Extrapolation  scheme 
crosses  from  zone2  into  zonel  during  the  10th  timestep  (starting  at  t  =  0.40  seconds). 

Table  [9]  shows  that  the  analytic  and  numerical  system  parameters  agree  to  approx¬ 
imately  15  decimal  places  at  the  end  of  timestep  9.  During  timestep  10,  each  simulation 
block  starts  in  zone  2  and  moves  to  zone  1.  It  is  important  to  note  that  all  sub-steps  of 
standard  Newmark  integration  conserve  energy,  including  those  that  cross  the  boundary 
between  zones.  Still,  the  extrapolation  of  displacements  and  velocities  fails,  and  is  ac¬ 
companied  by  a  significant  loss  in  accuracy  for  the  system  energy.  The  latter  results  are 
summarized  in  Table  [10]. 
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8  CONCLUSIONS 


This  paper  has  described  the  formulation  of  a  new  class  of  algorithms  for  integrating 
linear  second  order  equations,  and  those  containing  smooth  nonlinearities.  These  algo¬ 
rithms  are  based  on  ideas  taken  from  standard  Newmark  integration  methods,  and  extrap¬ 
olation  techniques.  For  the  algorithm  to  work,  the  underlying  Newmark  method  must  be 
stable,  second  order  accurate,  and  produce  asymptotic  error  expansions  for  the  response 
quantities  containing  only  even  ordered  terms.  It  has  been  shown  that  second  order  accu¬ 
racy  and  a  desirable  asymptotic  expansion  is  achieved  for  Newmark  parameters  7  =  1/2, 
and  any  value  of  f3.  This  means  that  (3  may  selected  on  the  basis  of  stability  considerations 
alone.  The  numerical  experiments  indicate  that  the  Newmark-Extrapolation  algorithm  is 
capable  (in  certain  situations)  of  tracking  response  quantities  in  excess  of  1,000,000  times 
more  accurately  than  the  standard  Newmark  method. 

The  current  plan  is  to  extend  this  work  in  two  directions.  First,  algorithms  of  the  type 
developed  here  lend  themselves  to  concurrent  computing.  A  prototype  environment  for 
distributed  numerical  computations  is  currently  being  developed.  The  algorithm’s  perfor¬ 
mance  in  this  environment  will  be  reported  in  future  papers.  Second,  it  is  planned  to  extend 
the  extrapolation  technique  to  single  rigid  body  and  coupled  rigid  body  applications.  The 
challenge  here  is  to  try  and  find  ways  of  doing  extrapolation  within  the  Lie  Algebra  of  the 
application  at  hand.  By  doing  this,  the  hope  is  that  it  may  be  possible  to  compute  refined 
estimates  of  system  response  without  destroying  invariants,  such  as  conservation  of  energy. 
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